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ABSTRACT 

The variation of radial impurity distribution induced by 
surface tension driven flow increases as the zone length decreases 
in silicon crystals grown by floating zone melting. In combined 
buoyancy driven and surface tension driven convection at the gra- 
vity of earth, the buoyancy contribution becomes relatively smal- 
ler as the zone diameter decreases and eventually convection is 
dominated by the surface tension driven flow (in the case of silicon, 
for zones of less than about 0. 8 cm in diameter). Preliminary 
calculations for sapphire suggest the presence of an oscillatory 
surface tension driven convection as a result of an unstable melt 
surface temperature that results when the zone is heated by a ra- 


diation heater. 
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SUMMARY 


This report is supplementary to the 1974 Annual Report [l] on the 
analysis of surface tension, driven convection in floating zone melting. 
The detailed mathematical development and numerical method used 
for calculation of the momentum, heat and mass transfer equations 
are reported in Reference [l], and are omitted in this report. 

The influence of the operating parameters such as zone geometry, 
heating methods (electron beam or radiation heating) and gravity of 
earth are investigated. We have attempted to solve transient equa- 
tions for floating zone melting of sapphire for a possible oscillatory 
convection with radiation heating. The following conclusions are made 
from the results of our computer calculations for silicon, except 
where noted; 

1) Radial distribution of impurity in the crystal increases 
as the zone length decreases due to the fact that there is 
less convection at the center of the short zone. 

2) As the zone size (length and/or diameter) decreases, the 
buoyancy effect due to the earth’s gravity diminishes -and 
convection is dominated by the surface tension induced 
flow. 

3) The floating zones of high Prantl number materials (for 
example, sapphire) are susceptible to an oscillatory con- 
vection due to the unstable temperature at the free melt 
surface when a radiation heater is employed. 
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4) The surface tension, driven convection tends to be more 
well defined with electron beam heating (i, e., concentra- 
ted heat supply) than with radiation heating (i, e, , distri- 
buted heat supply). 

5) The maximum stream function of pure surface tension 
driven flow increases as surface tension parameter M' 
increases with the relationship of ijr max = 0. 199 (M^)^* 

6) The buoyancy effect is relatively small near the free 
melt surface and large near the center of the zone for . 
combined surface tension driven and buoyancy driven 
flows. 

During the year, three papers were published [2-4] and a fourth 
was accepted for publication [5],. A paper was presented at the 
Stanford Conference on Crystal Growth in July [6] and a paper on 
related studies was presented at the Pasadena meeting of the Ameri- 
can Physical Society in December [7]. 
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I. INTRODUCTION 


As reported in the 1974 Annual Report [l], numerical analysis 
of convection in floating zone melting of silicon has shown that the 
surface tension driven convection is vigorous. 

During the current contract period (January 1975 to December 
1975), the influence of both operating parameters and property par- 
ameters on surface tension, driven convection were investigated. 

The following operative parameters were considered: 

1) Zone geometry, such as aspect ratio and zone diameter. 

2) Relative importance of buoyancy driven- convection due 
to the gravity of earth. 

3) Heating modes, including electron beam heating and ra- 
diation heating. 

In addition to silicon, sapphire was selected to investigate the 
influence of physical properties on the surface tension driven flow. 

On earth, the relative importance of the buoyancy effect is ex- 
pected to increase as zone size increases, since the volume of the 
melt increases at a greater rate than the free melt surface area. 
This argument was applied to a small zone, in which the surface 
tension driven flow is expected to be predominent due to the large 
free melt surface area relative to the melt volume. From the par- 
tial differential equation for the momentum transfer, a parameter 
combining the Grashof number (Gr^^) and the dimensionless surface 
tension parameter (M^ were derived in order to define the relative 
importance in the combined buoyancy driven and the surface tension 
driven flow. 
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The floating zones of high Prandtl number materials, such as 
organic compounds and metal oxides, are more susceptable to os- 
cillary convection when radiation heating is used, since the heat 
radiated on the free melt surface from the heater tends to accumu- 
late as the melt moves along the surface due to the low heat conduc- 
tion rate from the surfact into the hubs melt, A transient solution of 
the transport equations will be necessary to further investigate the 
possibility of surface tension driven oscillatory' convection in the 
floating zone melting of sapphire. 


II. NUMERICAL ACCURACY AND CONVERGENCE TEST 

The accuracy and convergence in the previous numerical solutions 
of the surface tension driven flow in. floating zone melting were exam- 
ined for a representative case as follows: 

1) The same number of grid points were used, but a higher num- 
ber of iterations was made by placing a more stringent accuracy limit. 
The results were compared with the previous results in order to verify 
the convergence, 

2) The grid size was reduced and the computation was carried 
out for a representative case with a sufficient number of iterations, 

A. Convergence Test 

Convergence tests were made by increasing the number of itera- 
tions until no appreciable difference was found in the solution of the 

test problem. Iteration was terminated when the difference between 

f s^l ) til 

the current value of stream function, . ’ at the (s+1) iteration 
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and the old value of stream function at the (s)^^ iteration passed 
the accuracy limit at every grid point. We tested the case shown in 
Figure 7 in the Fi-rst Annual Report [l] where the accuracy limit was 


J 


4 . 1 1 


< 0.005 


under which the required number of iterations was 217, When the 
accuracy limit of 0.001 instead of 0.005 was used, the required num- 
ber of iterations was 383, However, there was no essential differ- 
ence in the pattern of stream functions, the fact of which indicates a 
sufficient convergence in this problem. 

B. Grid Refinement Test 

The stream functions, temperature and concentration fields 
shown in the previous figures (Figures 7 and 18, Ref. [l] ) were 
tested as a representative case by refining the grid size from 11 x 21 
to 21 x41 grid points and comparing them. An accuracy limit of 
0.0005 was employed for the stream function. The computation of 
stream function was terminated after about 2100 iterations. The 
results in stream functions, temperature and concentration fields 
with 21 X 41 grid points are shown in Figures 1 through 3. In Table I, 
the stream functions, vorticity and fluid velocities are compared 
with the results computed with 11 X 21 grid points. The numerical 
value of the stream function showed about 1% deviation near the free 
melt surface where the convection is weakest. However, there was 
no essential difference in the pattern of the isotherms and the stream- 
lines, nor in the intensity of the flow, as shown in Figure 1. The 
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10.11 10.37 



10.70 



FREEZING INTERFACE 


Figure 2. 


Dimetisioaless impurity concentratiou fields, C=m/m^, 
in a silicou melt at steady state for a zone travel rate 
of = 5 cm/hr. with the same conditions and grid points 
as in Figure 1. k' =0. 1. 



AXIAL POSITION, 



DIMENSIONLESS IMPURITY CONCENTRATION 
{ C' = m/m| ) 


Figure 3. The calculated dimensionless impurity concentration 
profiles in the axial direction at R = 0, 0.5 and 1. 0 
under the same conditions and grid points as in Fig- 
ure 1. 
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location of the vortex center near the heat source was found to be in- 
corporated with the individual grid point nearest to the axial and the 
radial axis. 

The dimensionless concentration field of the impurity, C=m/m^, 
computed with 21 X41 grid points is shown in Figure 2 and is essen- 
tially the same as that in Figure 18 in the Annual Report [l], except 
for an increased smoothness in the isoconcentration lines. The nu- 
merical values of concentration deviate slightly in the lower cell and 
deviate about 20 to 25% in the upper cell. The axial dimensionless 
concentration profiles at R = 0, 0, 5, 1.0 are plotted in Figure 3. The 
wavy concentration profile at the melt free surface in previous Figure 
19 [l] becomes smoother, as shown in Figure 3. This indicates that 
the wavy profile was associated with small errors at the individual 
grid points due to the coarse grid size. Another possible cause of 
numerical discrepancy may come from the coherent dependence of the 
interfacial gradients on the grid size, because the gradients are rep- 
resented in a linear form in the finite difference equations. 

C. Discussion 

Grid refinement tests for the stream function, temperature and 
concentration fields for the case shown in previous Figures 7 and 18 
[l] were carried out by doubling the grid points from 11 X 21 to 21 X 41, , 
The required computer time with 21 X41 grid points was about 20 
times more than that with 11 X 21 grid points. Not only did no additional 
vortex appear in the grid refinement tests, but also no fundamental 
difference in the overall patterns was found, -although the concentration 
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profile and the vortex center near the heat sonrce tended to he asso- 
ciated with the individual grid points with the coarser grid system. 
Convergence was much faster with the coarse mesh size than the 
finer mesh size. 

We conclude that there is no significant advantage in using a 
finer mesh size at the expense of increased computer time, unless 
the flow pattern is very complex due to flow separation with many 
small vortices. 


m. INFLUENCE OF THE ZONE LENGTH ON IMPURITY 
SEGREGATION 


For an aspect ratio (X/a>) of other than one, the dimensionless 
surface tension parameter at the free melt surface, M^, should be 
defined in terms of the zone radius, a, and one half of the zone 
length, 1 , such as 


aP(T - T )a- 

“ TTZ 


U i 


Computations were performed for a silicon zone of 1 cm diam- 
eter with the following variables: 

1) A/a=-|, (T^- T^) = 0.0125°Q, and M' = 350. 

2) A/a = 2, {T^- T^) = 0.2°C, and M' = 350. 

The heating modes employed in the present calculations are: 

1) Parabolic surface temperature profile, which is assumed 
to be a radiant heating mode, 

2) Ring heat source at the center of the zone with electron 


beam heating 



The results for .the computed streamline, temperature and con- 
centration fields are shown in Figures 4 through 11, The vortex cells 
are stretched in the long zone (A /a = 2) and contracted in the short 
zone (jj/a = -g- ) and no additional multiple cells are generated. How- 
ever, in a short zone with electron beam heating, the secondary vor- 
tex cell, which appears in the long zone, is combined with the pri- 
mary vortex cell to make a single cell near the ring heat source, as 
shown in Figure 6, 

The impurity incorporation into the crystal is plotted in Figure 
12, It shows that the radial distribution of the impurity in the grown 
crystal is not sensitive to the mode of heating, .but is sensitive to the 
zone length. The radial inhomogeneity of impurity is greater for the 
crystal with the shorter zone, which is expected from the fact that the 
convection of the melt near the center of zone is relatively weak com- 
pared with that in a long zone. The isotherms near the melt/solid 
interfaces are curved in the short zone, while those in the long zone 
are nearly flat. This, in turn, reflects the tendency for the inter- 
faces to deviate from a flat to a curved shape as the zone length is de- 
creased. The characteristics of the flow patterns and the physical 
and operational parameters are summarized in Table I. 


IV. COMBINED SURFACE TENSION DRIVEN AND BUOYANCY 
DRIVEN CONVECTION 

As discussed earlier, is the primary factor characterizing 
the surface tension driven flow regardless of the size of the zone when 
the buoyancy effect is neglected. However, the buoyancy effect is ex- 
pected to increase as the diameter of the zone increases, since the 


13 




Figure 4. Computed dimensionless streamlines, ijf for surface 
tension driven flow in a floating zone at zero gravity 
with, a parabolic temperature profile on the free melt 
surface with (1^ - T^) = 0.0125°C, a=0,5 cm, 0, 25 
cm, M^-350, and •^ = 5 cm/hr. Calculations per- 
formed with 11 X 11 grid. 
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ure 5. Computed dimensionless impurity concentration fields, 
C =m/m^ at steady state in a silicon melt for the same 
conditions as in Figure 4 with Sc = 5. 0, = 0, 1, and 

D = 7. 6 X 10“ - cm®/sec. 
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Figure 6. Computer calculated isotherms, 0 (dotted lines) and 
streamlines, ^ (solid lines) for surface 'driven flow 
in floating zone melting at zero gravity, using the 
computed temperature fields for electron beam heat- 
ing of silicon with (T^ - Tj,)-= 0* 0125° C, Tc = Tn, 5g = 
0,3, a = 0,5cm, <1=0.25 cm, M^=350, Pr = 0.023, 

Vo = 5 cm/hr, and 0 = (T-T 5 j)/(To- TJ. Calculations 
performed with 11x11 grid. 
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Figure 7. Computed dimensionless impurity concentration fields, 
C =m/m^ at steady state in a silicon melt for the same 
conditions as in Figure 6 withSc = 5.0, k^ = 0.1, and 
D = 7. 6 X lO”^ cm2/sec. 
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Figure 9. Computed dimensioaless impurity concentration fields, C at steady state in a silicor 

melt forthe same conditions as in Figure 8, with Sc = 5. 0, k„= 0, 1, and D = 7. 6 x 10“‘^cm2/se< 





Figure 10, Computer calculated isotherms, 0 (dotted lines), and sti-eamlines, i|; (solid lines) for 

surface tension driven flow in floating zone melting at zero gravity, using the computed 
temperature fields for electron beam heating of silicon with (T^,- Tj,,) = 0. 2“C, Tc= T^, <?j= 
0.3, a = 0.5 cm, ^ = 1 cm, M'=350, Pr = 0.023, Vc=5cm/hr, and Q = (T-Tj,)/(Tg- ). 
Calculations performed with 11 X41 grid. 








DIMENSIONLESS IMPURITY CONCENTRATION IN 

CRYSTAL, C = m/m, 



DIMENSIONLESS RADIAL POSITION, R 


Figure 12. Radial impurity coacentration profiles in the crystal 
at steady state freezing rate of 5 cm/hr. calculated 
. from the interfacial concentration in Figures T, % 
and 11 for interfacial distribution coefficient of = 

0 . 1 . 
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buo/ancy driven, flow is due to the volume effect while the surface ten- 
sion driven flow is due to the free liquid surface area effect. ] 

In order to define a new parameter for combined surface tension 
driven and buoyancy driven flow, the dimensionless velocities of flow 
were redefined based on the reference velocity. The parameters de- 
rived from the momentum equation are and 1/M^ The com- 

bined surface tension driven and buoyancy driven flows are compared 

2 

with pure surface tension driven flow as a function of Gv-^/Mf and 
the geometry of zone. 


A. Model for Analysis of the Combined Convection 

In the present study, we derived parameters combining Gr^ and 

in order to obtain a better correlation for the combined surface 

tension driven and buoyancy driven flow. Employing the reference 

2 2 

velocity defined as v^ = ct(T^- T^)a /|aA , the surface tension par- 
ameter is absorbed in the momentum equation instead of in the 
free melt boundary condition. By defining t = tv /a, V = v /v , V = 

V /v . Z = z/a, and R = r/a, the following momentum equation is ob- 
z o 

tained through the similar operations made in the annual Progress 
Report: 


au) _ Gr^ ^ , S M 5(Rw) N't 

SR az 2 I 2 ^ ^R\R aR JJ 


where 


and 


-1 Sir _ 1 3i!f j,/_ 

r " R aZ" z'^RaR" . 2,2 


^^h = 


gP(T^-TJa 


2 


v‘ 
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2 

The parameter, , can be separated into the property parameter 

and the operation parameter as follows; 




« 2 2 
gpu P 
2 


} 


{ 


(T - T_)a' 
o m 


} 


property operating 

parameter parameter 


( 3 ) 


B. Influence of the Zone Length 

Inspecting Equation (3) the buoyancy effect is expected to increase 
as the zone length, H, increases. Calculations were performed for 
both short and long zones of silicon using the following values: 


Zone Diameter 
in cm 

Zone Length 
in cm 

{T - T ) 
^ o m^ 

in °C 

M 


2 

1 

0.00625 

350 

77.5 

2 

4 

0. 1 

350 

1240. 


The calculated streamlines and isotherms are shown in Figures 13 
and 14. For the identical value of = 350, the buoyancy effect was 
pronounced in a long zone and negligible in a short zone. 

C. Influence of the Zone Diameter 

2 

When the zone aspect ratio is one (i.e,, a = f), Gr^lMf be- 
comes proportional to a/(T^~ This suggests that the buouancy 

effect increases as the zone diameter increases and/or the tempera- 
ture difference in the zone decreases. Therefore, the relative-in- 
fluence of the buoyancy driven convection is expected to increase as 
the efficiency of thermal insulation of the zone increases. 
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Figure 13. Computer calculated isotherms, 0 (dotted lines) and 
streamlines, \]r (solid lines) for surface tension driv- 
en flow in floating zone melting at the gravity 'of earth, 
using the computed temperature fields for electron 
beam heating of silicon with (T^- T^) =0. 00625® C, Tg= 
Tjj, dj = 0.3, a = 1 cm, jJ=0.5cm, m'= 350, Pr = 0.023, 
Grj^ = 77.5, v^ = 5 cm/hr. , and 0 = (T-Tj,)/(To- TJ. 
Calculation performed with 21 X 21 grid. 
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Figure 14. Same conditions as in Figure 22 except (T^- Tj,) =0, 1“C, 
a = l cm, ^=2 cm, and Gr^ = 1240. Calculations per- 
formed with a 41 X 11 grid. 
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Both the combined convection and the pare surface tension 
driven convection were calculated with various zone diameters and 
and were compared in terms of the ratio of Nusselt numbers, 
r at the center of the freezing interface. That is 

1 , _ Nusselt No. at Ig _ 

Nusselt No. at Og ^^0 

The values used in the calculations for electron beam heated float- 
ing zone melting of silicon with the aspect ratio of one (a = jl) are 
shown in Table II. 

For silicon zones of 0. 2 cm to 8 cm diameter with M' = 350, 
earth's gravity and various ring heat source temperatures depending 
on the size of the zone, the convection patterns for combined buoy- 
ancy driven and surface tension driven flow as a function of zone di- 
ameter are shown in Figures 15 through 20. As the zone diameter 
increases, the vortex cell below the center of the zone gradually di- 
minishes while the convection of the vortex cell above is significantly ' 
enhanced. In order to relate the buoyancy effect to the zone diameter, 
the Nusselt number (i. e., dimensionless temperature gradient) as the 
solid/melt interface is calculated with and without the buoyancy ef- 
fect. At the center of the interface, the Nusselt number was most 
sensitive to the buoyancy effect and least sensitive at the periphery 
of the interface. The values of p with M^=350 are plotted as a func- 
tion of the zone diameter in Figure 21, Note that the buoyancy effect 
becomes negligible when the diameters of the zone are below about 
0 . 8 cm. 
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Table II. Ratio of Nusselt numbers at Ig and Og for 
various experimental conditions 


Zone Diameter 
in cm 

(T - T ) 
0 ra' 

in “C 

M' 


'^max 
at Og 

Nu /Nu- 
g 0 

= r 

2 

.0125 

175 

155 

4. 11 

. 3705/. 3689 

= 1.004 

4 

.00625 

175 

2480 

3.37 

. 7454/. 6623 

= 1. 125 

10 

.0025 

175 

3875 

4, 89 

1.220/. 8945 = 

^ 1.364 

0.2 

. 25 

350 

3. 1 

6. 64 

.3881/. 3880 

= 1.00 

0.5 

. 10 

350 

19.4 

6. 68 

.38487.3842 

= 1.002 

1 

.05 

350 

77.5 

6.84 

.3815/. 3765 

= 1.013 

2 

.025 

350 

310 

7. 19 

.3775/. 3521 

= 1.072 

4 

.0125 

350 

1240 

7.47 

. 4304/. 3626 

= 1. 187 

6 

.00833 

350 

2790 

7. 86 

. 5985/. 4551 

= 1.315 

8 

.00625 

350 

4960 

7.84 

. 7735/. 5626 

= 1. 375 

0. 2 

. 5 

700 

6. 2 

10. 76 

.3614/. 3604 

= 1.003 

0. 5 

, 2 

700 

38.8 

10. 22 

. 3648/. 3539 

= 1.031 

1 

. 1 

700 

155 

10. 10 

.4198/. 3402 

= 1.234 

2 

.05 

700 

620 

9.67 

. 5123/. 3282 

= 1.561 

4 

.025 

700 

2480 

11. 31 

. 5058/. 2982 

= 1.696 

6 

.0166 

700 

5580 • 

11.08 

. 5534/. 3094 

= 1.789 

. 0.5 

.30 


58. 1 

11. 64 

.3052/.-2966 

= 1.029 

1 

. 15 

1050 

232. 5 

11. 54 

. 3974/, 2936 

= 1.354 
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SOLID 


Computer calculated isotherms, 0 (dotted lines) and 
streamlines, ijf (solid lines) for surface tension driv- 
en flow in floating zone melting at the gravity of earth, 
using the computed temperature fields for electron 
beam heating of silicon with (T^- TJ =0. 25“ C, T<.= X,, 

= 0.3, a=;J=0.1cm, M'=350, Pr = 0.023, Gr^ = 3. 1, 
v« = 5cm/hr, and 0 = (T-TJ/(T,-* TJ. 





VACUUM 


Figure 17. 


Same conditions as in Figure 15 except (Tg- T ) = 0. 1°C, 
a =i, = 1 cm and Gr^= 310. " 
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Figure 18, Same conditions as in Figure 15 except (T^,- - 

0,0125“C, a = 2 cm and Grjj = 1240. 
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Figure 19. Same conditions as in Figure 15 except fT - T ^ 

0.00833“C, a = i =3 cm and Gr^= 2790. “ “ 










Figure 21. The ratio of the NusseXt number of combined surface tension driven and buoyancy driven 
flow at 1 g to that of pure surface tension driven flow at the center of the interfaces ver- 
sus the zone diameter for the floating zone melting of silicon with electron bearn heating, 
To=Tj,, c?g=0,3, M^=350, Pr = 0.023, and Vc= ,6 cm/hr. 


u> 

Ln 




Figures 2Z and 23 illustrate the effect of- increased zone diameter 
on the pure surface tension deiven flow at zero g for a silicon zone of ’ 

4 cm diameter with i = 2 cm, (T - T )- = 0.0125°'C and. M^=350, Both 
a parabolic surface temperature profile and a- ring heat source at the 
center of zone were employed. The computed streamlines with para- 
bolic temperature profile- are shown in Figure 22, The simultaneous 
solution for streamlines and temperature fields with electron beam 
heater are shown in Figure 23, The convection patterns shown in 
Figures 22 and 23 are similar to those obtained with a smaller zone 
diameter, shown in Figures 3 and 7 of the First Annual Progress 
Report [l]. The similarity in convection patterns- indicates that 
is the primary factor characterizing the surface tension driven flow 
regardless of the size of the zone for a given material. It should be 
noted, however, that as the diameter of the zone increases, will 
actually increase in proportion to about the square of the increase in 
the zone diameter, since the temperature difference- (T^- T^) will 
also increase in proportion to the increase in zone diameter. There- 
fore, the increase in zone diameter from I to 4 cm has the same ef- 
fect of increasing from 350 to about 1400, and convection would 
become stronger in the larger zone. 

The maximum stream functions for pure surface tension driven 

flow, ilr , are plotted as a function of in Figure 24. The rela- 
■max 

tionship is correlated by the' following expression: 


(M^) 


0.6 


( 5 ) 
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VACUUM 



Figure Z2, Computed dimensionless streamlines, ilr, for surface 
tension driven flow in floating zone melting of silicon 
at zero gravity, with a parabolic temperature profile 
on the free melt surface and (T,~ Tjj) =0, OIZS^C, a = 

A = 2 cm, 350, and v,. = 5 cm/hr. Calculations 
performed with 11 X21 grid. 
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Figure 23. 


Computer calcalated isotherms, 0, (dotted lines) and 
streamlines, \[f, (solid lines) for surface tension 'driv- 
en flow in floating zone melting at zero gravity, using 


the computed temperature fields for electron beam 
heating of silicon with (T^-TJ =0.0125“C, T„= T , S^- 
0,3, a=A = 2cm, M^=350, Pr = 0.023, Vc=5cm/hr., 
and 0 = (T-Tg) /(Tj,- T^). Calculations performed with 
11x21 grid. 


max 



Figure 24. Correlation of the surface tension parameter, M^, 
and the maximum stream function, of the sur- 

face tension driven flow in floating zone melting of 
silicon at zero gravity. The diameters of the zone 
in cm are: D = 0. 2(0), 1(©), 2(a), 4(A)< 8(x), 10(:5;'). 
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The values of T are plotted as a function of Gr^^/M in 

Figure 25. As the value of Gr^/M'^ decreases, the buoyancy effect 

diminishes and the value of LogP approaches zero, 1/M^ and 
2 

are shown to be independent parameters defining the con- 
vection. This is also indicated by the momentum Equation (2) which 

2 

contains both parameters, and 1/M^ No significant im- 

provement was found when the values of p were plotted as a function 
of Grj^/M' instead of Gv^Im!^, 

Finally, the buoyancy effect is represented in terms of the flow 
components in the zone in Figure 26. The axial flow velocity increases 
in the direction of gravity as the Grashof number increases. When 
the Grashof number is 6.2, the flow patterns approach those for pure 
surface tension driven flows. 

V. FLOATING ZONE MELTING OF SAPPHIRE AT ZERO G 

From the survey of the properties of materials of possible in- 
terest for space processing, sapphire was selected as the candidate 
material in view of its commercial importance and its unique transport 
properties, which are considerably different from silicon. The phys- 
ical data surveyed for Al^O^ a-re substantially different from those es- 
timated in the Annual Report [l]. The new data for listed in 

this report, were used for the present study. . 

A. Computations with Free Parameters 

The streamlines and the temperature profiles in the zone for 
electron beam heating were computed with free parameters chosen 
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LOG (T) 



LOG (Gr^/M-) 


Figure 25, The ratio of Nusselt numbers, P defined by the equation (4), for the combined surface 
tension driven and buoyancy driven flow at 1 g to that for pure surface tension driven 
flow at zero g at the center of the freezing interface (bottom interface) versus Gru/M'®. 





AXIAL POSITION, 


TOP INTERFACE 


= 2480 ^ 



BOTTOM INTERFACE 
I I 


Figure 26. The influence of the Grashof number, Grj^ in the 
combined surface tension driven and buoyancy- 
driven flo-w in the floating zone melting of silicon 
on the axial velocity- component of flow, along 
the center of zone for M^=700, {- indicates the 

same direction as the gravity, and -f- the opposite) 



from properties intermediate between those of silicon and sapphire, 

> 

such as, ^ = 0. 7, k = 0.024 cal/cm*C, C = 0, 3 cal/g°C, p^ = 

P * 

.3. 05 g/cm^, p, = 2, 72 CentiPoise, ^^ = 0.06 djrne/cm®C, Pr = 0.34, 

(T - T ) = 0.05'‘C, and = 22, The results are shown in Figure 27. 
' o rn 

From the results studied so far with variable properties of 
materials we conclude that the intensity of the surface tension dri- 
ven flow is primarily dependent upon the surface tension parameter 
M', and the flow ,-patter ns change mainly as a result of variations in 
the Prandtl number. 

B, Survey of the Physical Properties of Sapphire 

Most of the sapphire melt physical properties required for the 
present study are available in the literature. The data employed are 
listed below: 

1) Temperature coefficient of surface tension, By/ST, 0.465 
dyne/cm°C [8]. 

3 

2) Density, p, 3,05 g/cm 

3 

3,66 g/cm (solid). 

3) Viscosity, p, 0.6 g/cm sec [8], 

4) Specific heat, C^, 0.3 cal/g°C [9]. • 

5) Thermal conductivity, k, 0.015 cal/cm°C sec [9]. 

6) Emissivity, S > 0.5 [9]. 

Several potential problems are encountered in the floating zone 
melting of sapphire, particularly, 

1) Suitability of electron beam heating, and 

2) Effect of internal thermal radiation. 
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Figure 27. Computer calculated Lsotheims, 0, (dotted lines) 
and streamlines, ijf, (solid lines) for surface ten- 
sion driven flow in floating zone melting at zero 
gravity, using the computed temperature fields 
for electron beam heating with free parameters 
such as: (T^- TJ =0.05'’C, T<,= T„, (S = 0.7, a=A 
= 0,5 cm, k = 0.024 cal/cm®C sec, p^=3,.05 g/cm®, 
ji = 2. 72 centipoise, M'= 22, Pr = 0,34, v<;= 5 cm/hr, 
and 0 = (T-Ta)/ (T^- TJ. Calculations performed 
with 11x21 grid. 




-2 

The electrical conductivity drastically increases from .10' 
to 15 ohm ^ cm ^ when AI 2 O 2 is melted [lO], Therefore, electron 
beam heating maybe feasible if the material is preheated to the 
melting point by other means. The infrared absorption of an AI 2 O 2 
melt is approximately 12 times that of the solid at melting point 
[11.], which indicates that the radiation effect is relatively less im- 
portant in the melt than in the solid. At present, it does not seem 
practical to rigorously include radiation in heat transfer calcula- 
tions, due to the complexity of the problem. Therefore, approxi- 
mate computations were made by using the apparent thermal con- 
ductivity, which takes the radiation heat transfer effect into account 
in terms of the conduction heat. Two modes of heating are consid- 
■ered in the following. 

C. Computations with Electron Beam Heating 

In order to determine the required heat shield temperature^ 

T^, for a given temperature of ring heat source, T^, the tempera- 
ture fields in the ingot were computed for various heat shield tem- 
peratures. For (Tq- T^) = 180® C and 72° C, the required values of 
T^ to produce a zone of length equivalent to diameter were approxi- 
mately E033°C and 2010° C, respectively. The temperature profiles 
in the ingot are shown in Figures 28 and 29. From the physical data 
for sapphire, was estimated to be 35 with = 18°C, and 

140 with (T - T ) = 72° C. The estimated Prandtl number is 12, 
o m 

and the radiation shape factor ^ is 0.25 (see the following section). 
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In the zone of a=j6 =0.5 cm the computed streamlines and tempera- 
ture fields are shown in Figure 30 for M' = 35 and in Figure. 31 for 
= 140, When = 140, the vortex near the ring heat source was 
larger compared with that in the silicon zone, 

D, Computations with Radiation Heating 

In order to simplif/ the problems of radiation heat transfer 
between the heater and the material, the following assumptions were 
made; 

1) The heater temperature is uniform at 

2) The radius of the heater approaches that of the zone. 

The radiation shape factor, ^ , between infinite parallel plates can 
be used under the condition of 2). For example, the heat flux be- 
tween the heater and the zone, q, is 

q = a^(T^- T^) (6) 

where — 1 and S and are the emissivity of the ma- 

< ' O O u S il 

S n 

terial in the zone and that of the heater, respectively, and a is the 
Stefan-Boltzmann constant. 

The emissivity [9] of metals of high melting point (for example, 

W, Mo, Ni, and Ta) is about 0. 33 at 2600 2800 ®C which is close to 

the required heater temperature for the floating zone melting of sap- 
phire, With 0. 33 and 3^ = 0. 5, the estimated J from Equation 
(6) is 0.25. 
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Figure 30. Computer calculated isotherms, Q,. {dotted lines) and 
streamlines, ijr, (solid lines) for surface tension dri- 
ven flow in floating zone melting at zero gravity, us- 
ing the computed temperature fields for electron beam 
heating of sapphire, with {T^- TJ = IS** C, T.= 2037“ C, 
To=2033“C, J'^O.ZS, k = 0.015 cal/cm“C sec, a=^ = 
0.5 cm, M^=35, Pr = 12, v<.= 5 cm/hr., and 0' = (T-TJ/ 
(T^- T|j). Calculations performed with 11 X 21 grid. 



Figure 31, Computer calculated isotherms, 0, (dotted lines) and 
streamlines, i|f, {-solid lines) for surface tension driv- 
en flow in floating zone melting at zero gravity, using 
the computed temperature fields for electron beam 
heating of sapphire, with (T,- T„) = 72* C, T3=2037'’C, 
To= 2010*C, ^=0,25, k = 0, 0 15 cal/cm“ C sec, a=^ = 
0.5 cm, M'= 140, Pr=12, v<.= 5 cm/hr., and 0 = (T-Tg) / 
(T - T ). Calculations performed with 11 X 21 grid. 



The temperature fields in the sapphire ingot without convection 
were computed, for radiation heating and are plotted in Figure 32, 

The required heater temperature without convection was 2637* C. 
With the convection in the zone, however, the heater temperature 
can be lowered to 2337°C because of the more effective heat trans- 
fer to the convective zones. 

With the radiation heater applied to sapphire floating zone 
melting, we were not able to obtain a steady state solution, since 
the temperature of the free melt surface did not reach a stable value. 
The temperature gradient near the center is reversed slightly with 
the increased number of iterations, which indicates that the flow 
might become oscillatory. Therefore, we feel that it is of interest 
to investigate the transient solutions for the floating zone melting of 
sapphire with a radiation heater. 


E. Oscillatory Convection with Radiation Heating 

Transient solutions are required to investigate the oscillatory 
flows in the floating zone melting of sapphire when thermal radiation 
heating is employed. The equations of momentum and heat transfer 
are developed for transient solutions as follows; 

1) Equations for Transient Convection 

The dimensionless momentum equation for pure surface 
tension driven flow is, 


, 

bT 


b(V^Uj) 

" 11 “ 


, _b_ f ib(Rtu) \ 

bZ bH RbS. J 

bZ 


( 7 ) 


The dimensionless heat transfer equation is. 
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where, t = — ^ . V = * V = — , and other param- 

’ ' Z ' r V z V 
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eters are defined as previously. 

The approximate finite difference equations employed 
for the time dependent terms, 5 u)/5t and B0/5T are; 


<u 



(9) 


M = . (10) 

5t At 
2 

Where At = vAt/a and the terms with the asterisk are 
the values at the end of time stem. At* The same nu- 
merical schemes as used in the steady state calcula- 
tions are employed for the other terms in the equations. 


2) Numerical Computations 

The initial values for the stream function and the tem- 
perature were taken from one of the quasi steady state 
solutions. We were not able to obtain a converged solu- 
tion with the time increment, both in the temperature 
and in the stream functions. 


When the time dependent term was excluded from the momentum 
transfer equation (i. e., the transient condition is included only in the 
heat transfer calculation), the computer solutions were converged 
with a result that showed the oscillatory convection nature. 
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The calculation procedure used is as follows: 

1) Initialization of the temperature field from the quasi 
steady state solution. 

2) The stream function (based on steady state analysis) 
was computed using the initial temperature field from 
step 1). 

3) The temperature field at the end of the first time in- 
crement of was computed using the flow field from 
step 2). 

4) The stream function (based on steady state analysis) 
was computed using the temperature field from step 
2 ). 

5) Steps 3) and 4) were repeated. 

As an illustration, the stream functions and temperature 
fields were calculated in the radiation heated floating zone melting 
of sapphire with the previously mentioned physical configurations 
and the time increment of At = 0*0393 (of ,it = 0.05 sec for" sajpphire). 
The results are shown in Figures 33_ through 37 in sequence of the 
time interval of 2^x a period of approximately one oscillation 
cycle. The corresponding temperatures at the free melt surface 
are depicted in Figure 38. The convection cycle shown above was 
repeated with an extended time. 
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SOLID 


Computer calculated isotherms, 0, (dotted lines) and 
streamlines, ilr, (solid lines) for surface tension dri- 
ven flow in floating zone melting at zero gravity, us- 
ing computed temperature fields at t = 0 for radiation 
heating of sapphire, with (T^- T,,) = 300°C, T5j=2037°C 
To=23“C, k=0.015 cal/cm°C sec, J-=0.25, a=^=0.5 
cm', M'=585, Pr = 12, Y:= 5 cm/hr, and0=(T-Ta)/ 

(T,j- Tg). Calculation performed with 11 X 21 grid. 













AXIAL POSITION, 



DIMENSIONLESS TEMPERATURE, 0 = (T-Tm)/(Th-Tm) 


Figure 38, Temperature oscillation at the free melt surface. The numbers 1 through 5 indicate 
the temperature profiles corresponding to Figure 33 through 37 , in sequence. 
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VI. PROPOSED ADDITIONAL STUDIES 

The silicon calculations have demonstrated that surface tension 
driven convection can be significant even at the gravity of earth. Pre- 
liminary calculations on sapphire indicate' that convection may be os- 
cillatory when the molten zone is maintained by radiation heating. It 
is desirable to expand the studies on sapphire to further examine the 
question of oscillatory convection and to extend the analysis of sur- 
face driven convection to other systems of interest to space proces- 
sing. Calculations of the type recently completed on InSb [12] and in 
progress on InSb-GaSb alloys [6, 13] should be expanded and extended 
to other solid solution systems of potential commercial value, 

A number of experimental and theoretical studies have been 
conducted on various aspects of convection as related to space pro- 
cessing, The degree of interaction between investigators has not 
been optimum in all cases and all areas of conflict between theory 
and experiment have not been resolved. This would appear to be an 
excellent time to critically review the work that has been completed 
to date and determine the nature of future studies that should be con- 
ducted both to resolve conflicts and to maximize the benefits of pre- 
vious investigations. 
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APPENDIX A 


DEFINITION OF SYMBOLS 


a Radius of the molten zone, cm 


b 



C 


C 


ct 



D 

g 

§e 


Gr 


m 




k 



SL 


Distance of the top surface of the solid from the center 
of the zone, cm 

Specific heat, cal/g®K 

Dimensionless impurity concentration in the liquid, 
defined by m/m^ 

Dimensionless impurity concentration in the feed 
crystal, defined by m^^/m^ 

Dimensionless impurity concentration in the liquid at 
the crystallizing interface, defined by m^^ b^^t 

Dimensionless impurity concentration in the liquid at 
the melting interface, defined by m. 

2 

Binarymolecular diffusion coefficient, cm /sec 

,2 

Acceleration, cm/sec 

2 

Earth's gravitational field, 980 cm/sec 


Grashof number for heat transfer (/x ratio of the tempera- 
ture difference buoyancy force to the viscous force) de- 
fined by g3(T^- 

Grashof number of mass transfer (~ ratio of the composi- 
tion byoyancy force to the viscous force) defined by 
gam^a'^/v^ 

Integer denoting grid station in axial direction, i = l,2, 3, ... 
Integer denoting grid station in radial direction, j = 1, 2, 3, . . . 
Thermal conductivity, cal /cm sec®C 


Interfacial distribution coefficient (ratio of the impurity 
concentration in the crystal to that in the liquid at the 
interface) 

One half of the liquid zone length, cm 
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m 


m. 


m 

m 

m 


cb 

ct 

i,b 


m 




3 

Impurity concentration in liquid, b/cm 

Impurity concentration the feed solid would have if it 
were melted, g/cm^ 

3 

Impurity concentration in the crystal, g/cm 

3 

Impurity concentration in the feed crystal, g/cm 

Impurity concentration in the liquid at the crystallizing 
interface, g/cm^ ' 

impurity concentration in the liquid at the melting in- 
terface, g/cm3 


M' 


Nu 



Dimensionless surface tension parameter at the free 
liquid surface (~ ratio of surface tension force to vis- 
cous force) defined by pct(T - T 

V r \ o rn. ^ 

Nusselt number (dimensionless temperature gradient, 
or ratio of total heat transport to heat condiction) de- 
fined by |90/hz|±jj/a 

Nusselt number at the center of interface at the gravity 
of earth 


Nuq Nusselt number at the center of interface at zero 
gravity 

z 

p Dynamic pressure in fluid, dyne/cm , P-p 

2 

P^ Static pressure in fluid, dyne /cm , defined by p^g(;J-z) 

2 

P Total pressure in fluid, dyne/cm 

Pr Prandtl number,|j,Cp/k (ratio of the momentum diffusivity 

to the thermal diffusivity) 

r Radial coordinate, measured from the center of the zone, 

cm 


R Dimensionless radial coordinate, measured from the center 

of the zone, r/a 

AR Dimensionless radial distance between grid points 

Sc Schmidt number (ratio of the momentum diffusivity to the 

molecular diffusion coefficient), ^i/pD 
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Sh 


Sherwood number fdimensionless concentration gradient), 
defined b 7 j3C/9Z|±jS/a 

t Time, sec 

T Temperature, °K 

T Temperature of the middle of liquid zone at the periphery, 

° “K 

T^ Temperature of the surroundings, ®K 

T^ Interfacial temperature (melting point), 

Tj^ Temperature of the heater, ®K 

T Either (T^- T^) if 0 is defined by the Equation, 

9 = (T-T^)/{T^- T^), (T^- T^) if 0 is defined by the 

Equation, 0 = (T-T^) /(T^- T^), or ( ^ T^) if© is 

defined by the Equation, 0 = (T-T )/(T,“ T ) 

m' n m 

Zone travel rate, cm/sec 

2 2 

Reference velocity, defined by cr(T^- 

V Axial velocity component of the fluid, cm/sec 
z 

V, Dimensionless zone travel rate, v^a/v 

V Dimensionless radial velocity component of the fluid, 

^ v^a/v _ 

Dimensionless ctxial velocity component of the fluid, 
v^a/v 

Dimensionless interfacial flow rate, ^,,{P<,/pf) 

z Axial coordinate, measured fromdie center of the zone, 

cm 

Z Dimensionless axial coordinate, measured from the 

center' of the zone, z/a 

AZ Dimensionless axial distance between grid points 

a (Bp/Sm)„ -o/p] volume expansion coefficient due to a con- 

centration change in fluid, l/(g/cm'^) 

3 0p/3T) _/pJ volume expansion coefficient due to a tem- 

m, p 

perature change in fluid, 1/°C 
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y 

T 


j 

0 


V 

P 

Pc 

pf 

G 

a 

T 

ilr 




(JU 


Surface tension between melt and vapor, d 7 ne/cm 

Ratio of tbe Nusselt numbers for combined surface ten- 
sion driven and buoyancy driven convection, defined by 

Emissivity of thermal radiation 
Geometry factor for radiation heat transfer 
Dimensionless temperature, (T-T^)/AT 

Dimensionless temperature of the surroundings, 

(T -T* )/AT 
c m 

Viscosity of the fluid, g/cm sec 
Kinematic viscosity, (j./p, cm^/sec 

3 

Density of the fluid, g/cm 

3 

Density of the crystal, g/cm 

3 

Density of the fluid at melting point, g/cm 

-5 2 

Stefan- Boltzmann constant, 5. 668 X 10 erg/cm sec 
= k4 


Temperature coefficient of the surface tension, 
dynes/cm®K 


Dimensionless time defined by vt/a 


2 


Dimensionless stream function, defined by V = 




IStlf 

R3Z 


Dimensionless maximum stream function in the zone for 
pure surface tension driven flow 

Dimensionless vorticity, defined by yj = ^ 

1 M- + \ 
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